home *** CD-ROM | disk | FTP | other *** search
/ Acorn RISC PD-CD 1 / Acorn RISC PD-CD 1.iso / utilities / _graphics / graphics / _jpeg / c / jfwddct < prev    next >
Encoding:
Text File  |  1991-10-28  |  5.0 KB  |  180 lines
  1. /*
  2.  * jfwddct.c
  3.  *
  4.  * Copyright (C) 1991, Thomas G. Lane.
  5.  * This file is part of the Independent JPEG Group's software.
  6.  * For conditions of distribution and use, see the accompanying README file.
  7.  *
  8.  * This file contains the basic DCT (Discrete Cosine Transform)
  9.  * transformation subroutine.
  10.  *
  11.  * This implementation is based on Appendix A.2 of the book
  12.  * "Discrete Cosine Transform---Algorithms, Advantages, Applications"
  13.  * by K.R. Rao and P. Yip  (Academic Press, Inc, London, 1990).
  14.  * It uses scaled fixed-point arithmetic instead of floating point.
  15.  */
  16.  
  17. #include "jinclude.h"
  18.  
  19.  
  20. /* The poop on this scaling stuff is as follows:
  21.  *
  22.  * Most of the numbers (after multiplication by the constants) are
  23.  * (logically) shifted left by LG2_DCT_SCALE. This is undone by UNFIXH
  24.  * before assignment to the output array. Note that we want an additional
  25.  * division by 2 on the output (required by the equations).
  26.  *
  27.  * If right shifts are unsigned, then there is a potential problem.
  28.  * However, shifting right by 16 and then assigning to a short
  29.  * (assuming short = 16 bits) will keep the sign right!!
  30.  *
  31.  * For other shifts,
  32.  *
  33.  *     ((x + (1 << 30)) >> shft) - (1 << (30 - shft))
  34.  *
  35.  * gives a nice right shift with sign (assuming no overflow). However, all the
  36.  * scaling is such that this isn't a problem. (Is this true?)
  37.  */
  38.  
  39.  
  40. #define ONE 1L                  /* remove L if long > 32 bits */
  41.  
  42. #ifdef RIGHT_SHIFT_IS_UNSIGNED
  43. #define LG2_DCT_SCALE 15
  44. #define RIGHT_SHIFT(_x,_shft)   ((((_x) + (ONE << 30)) >> (_shft)) - (ONE << (30 - (_shft))))
  45. #else
  46. #define LG2_DCT_SCALE 16
  47. #define RIGHT_SHIFT(_x,_shft)   ((_x) >> (_shft))
  48. #endif
  49.  
  50. #define DCT_SCALE (ONE << LG2_DCT_SCALE)
  51.  
  52. #define LG2_OVERSCALE 2
  53. #define OVERSCALE (ONE << LG2_OVERSCALE)
  54.  
  55. #define FIX(x)  ((INT32) ((x) * DCT_SCALE + 0.5))
  56. #define FIXO(x)  ((INT32) ((x) * DCT_SCALE / OVERSCALE + 0.5))
  57. #define UNFIX(x)   RIGHT_SHIFT((x) + (ONE << (LG2_DCT_SCALE-1)), LG2_DCT_SCALE)
  58. #define UNFIXH(x)  RIGHT_SHIFT((x) + (ONE << LG2_DCT_SCALE), LG2_DCT_SCALE+1)
  59. #define UNFIXO(x)  RIGHT_SHIFT((x) + (ONE << (LG2_DCT_SCALE-1-LG2_OVERSCALE)), LG2_DCT_SCALE-LG2_OVERSCALE)
  60. #define OVERSH(x)   ((x) << LG2_OVERSCALE)
  61.  
  62. #define SIN_1_4 FIX(0.7071067811856476)
  63. #define COS_1_4 SIN_1_4
  64.  
  65. #define SIN_1_8 FIX(0.3826834323650898)
  66. #define COS_1_8 FIX(0.9238795325112870)
  67. #define SIN_3_8 COS_1_8
  68. #define COS_3_8 SIN_1_8
  69.  
  70. #define SIN_1_16 FIX(0.1950903220161282)
  71. #define COS_1_16 FIX(0.9807852804032300)
  72. #define SIN_7_16 COS_1_16
  73. #define COS_7_16 SIN_1_16
  74.  
  75. #define SIN_3_16 FIX(0.5555702330196022)
  76. #define COS_3_16 FIX(0.8314696123025450)
  77. #define SIN_5_16 COS_3_16
  78. #define COS_5_16 SIN_3_16
  79.  
  80. #define OSIN_1_4 FIXO(0.707106781185647)
  81. #define OCOS_1_4 OSIN_1_4
  82.  
  83. #define OSIN_1_8 FIXO(0.3826834323650898)
  84. #define OCOS_1_8 FIXO(0.9238795325112870)
  85. #define OSIN_3_8 OCOS_1_8
  86. #define OCOS_3_8 OSIN_1_8
  87.  
  88. #define OSIN_1_16 FIXO(0.1950903220161282)
  89. #define OCOS_1_16 FIXO(0.9807852804032300)
  90. #define OSIN_7_16 OCOS_1_16
  91. #define OCOS_7_16 OSIN_1_16
  92.  
  93. #define OSIN_3_16 FIXO(0.5555702330196022)
  94. #define OCOS_3_16 FIXO(0.8314696123025450)
  95. #define OSIN_5_16 OCOS_3_16
  96. #define OCOS_5_16 OSIN_3_16
  97.  
  98.  
  99. INLINE
  100. LOCAL void
  101. fast_dct_8 (DCTELEM *in, int stride)
  102. {
  103.   /* tmp1x are new values of tmpx -- flashy register colourers
  104.    * should be able to do this lot very well
  105.    */
  106.   INT32 tmp0, tmp1, tmp2, tmp3, tmp4, tmp5, tmp6, tmp7;
  107.   INT32 tmp10, tmp11, tmp12, tmp13, tmp14, tmp15, tmp16, tmp17;
  108.   INT32 tmp25, tmp26;
  109.   INT32 in0, in1, in2, in3, in4, in5, in6, in7;
  110.   
  111.   in0 = in[       0];
  112.   in1 = in[stride  ];
  113.   in2 = in[stride*2];
  114.   in3 = in[stride*3];
  115.   in4 = in[stride*4];
  116.   in5 = in[stride*5];
  117.   in6 = in[stride*6];
  118.   in7 = in[stride*7];
  119.   
  120.   tmp0 = in7 + in0;
  121.   tmp1 = in6 + in1;
  122.   tmp2 = in5 + in2;
  123.   tmp3 = in4 + in3;
  124.   tmp4 = in3 - in4;
  125.   tmp5 = in2 - in5;
  126.   tmp6 = in1 - in6;
  127.   tmp7 = in0 - in7;
  128.   
  129.   tmp10 = tmp3 + tmp0 ;
  130.   tmp11 = tmp2 + tmp1 ;
  131.   tmp12 = tmp1 - tmp2 ;
  132.   tmp13 = tmp0 - tmp3 ;
  133.   
  134.   /* Now using tmp10, tmp11, tmp12, tmp13 */
  135.   
  136.   in[       0] = UNFIXH((tmp10 + tmp11) * SIN_1_4);
  137.   in[stride*4] = UNFIXH((tmp10 - tmp11) * COS_1_4);
  138.   
  139.   in[stride*2] = UNFIXH(tmp13*COS_1_8 + tmp12*SIN_1_8);
  140.   in[stride*6] = UNFIXH(tmp13*SIN_1_8 - tmp12*COS_1_8);
  141.   
  142.   tmp16 = UNFIXO((tmp6 + tmp5) * SIN_1_4);
  143.   tmp15 = UNFIXO((tmp6 - tmp5) * COS_1_4);
  144.   
  145.   /* Now using tmp10, tmp11, tmp13, tmp14, tmp15, tmp16 */
  146.   
  147.   tmp14 = OVERSH(tmp4) + tmp15;
  148.   tmp25 = OVERSH(tmp4) - tmp15;
  149.   tmp26 = OVERSH(tmp7) - tmp16;
  150.   tmp17 = OVERSH(tmp7) + tmp16;
  151.   
  152.   /* These are now overscaled by OVERSCALE */
  153.   
  154.   /* tmp10, tmp11, tmp12, tmp13, tmp14, tmp25, tmp26, tmp17 */
  155.   
  156.   in[stride  ] = UNFIXH(tmp17*OCOS_1_16 + tmp14*OSIN_1_16);
  157.   in[stride*7] = UNFIXH(tmp17*OCOS_7_16 - tmp14*OSIN_7_16);
  158.   in[stride*5] = UNFIXH(tmp26*OCOS_5_16 + tmp25*OSIN_5_16);
  159.   in[stride*3] = UNFIXH(tmp26*OCOS_3_16 - tmp25*OSIN_3_16);
  160. }
  161.  
  162.  
  163. /*
  164.  * Perform the forward DCT on one block of samples.
  165.  *
  166.  * Note that this code is specialized to the case DCTSIZE = 8.
  167.  */
  168.  
  169. GLOBAL void
  170. j_fwd_dct (DCTBLOCK data)
  171. {
  172.   int i;
  173.   
  174.   for (i = 0; i < DCTSIZE; i++)
  175.     fast_dct_8(data+i*DCTSIZE, 1);
  176.  
  177.   for (i = 0; i < DCTSIZE; i++)
  178.     fast_dct_8(data+i, DCTSIZE);
  179. }
  180.